#include "TFile.h"
#include "TTree.h"
#include "TRandom.h"
#include "TMath.h"
#include "TProfile.h"
#include "TLorentzVector.h"
#include <iostream>
#include "TCanvas.h"
#include "TChain.h"
#include "TH2.h"
#include "TLorentzVector.h"
#include <iostream>
#include <complex>


using namespace std;

#include <gsl/gsl_sf_coulomb.h>
#include <gsl/gsl_sf_hyperg.h>
// double gsl_sf_hyperg_1F1(double a, double b, double x)
// int gsl_sf_hyperg_1F1_e(double a, double b, double x, gsl_sf_result *result)

int Coulomb_FG(double l,double eta,double x,double &F,double &G){
  gsl_sf_result F_gsl,Fp, G_gsl, Gp;
  double exp_F,exp_G;
  int status = gsl_sf_coulomb_wave_FG_e(eta,x,l,0,&F_gsl,&Fp,&G_gsl,&Gp,&exp_F,&exp_G);
  F = F_gsl.val;
  G = G_gsl.val;
  return status;
}

// double kr;
// double eta = alpha*m/(hbar*k);
// double F,G;
// Coulomb_FG(0,eta,kr,F,G);


// // double m1 = 1115.68;    //MeV
// double m1 = 938.27;    //MeV
// double m2 = 938.27;     //MeV
double m1 = 0.938272;    //GeV
double m2 = 0.938272;     //GeV
double m = m1*m2/(m1+m2);
int Z = 1;//Z1*Z2
// double hbar = 197.3269631; // MeV*fm/c
double hbar = 0.1973269631; // GeV*fm/c

double alpha = 0.0072973525693; // 1/137.035999084
double Euler_C = 0.57721566490153286;
double TwoPi = TMath::TwoPi();


double h_function(double eta){
  double h_value = -Euler_C-log(eta);
  double dh = 1;
  for(int n=1;dh<0.000001;n++){
    double dh = eta*eta/(n*(n*n+eta*eta));
    h_value = h_value + dh;
  }
  return h_value;
}

double norm2(complex<double> f){
    double n = f.real()*f.real() + f.imag()*f.imag();
    // return sqrt(n);
    return n;
}

complex<double> inverse(complex<double> a){
    return complex<double>(a.real(),-a.imag())/norm2(a);
    // return complex<double>(1,0)/a;
}

void complex_print(complex<double> a){
    cout<<Form("(%.20f,%.20f)\n",a.real(),a.imag());
}
void complex_print(complex<long double> a){
    cout<<Form("(%.20Lf,%.20Lf)\n",a.real(),a.imag());
}

///// 适用于|z|<30
complex<double> F_confluenthyper(complex<double> a,complex<double> z){
    complex<double> f(1,0);
    complex<double> df(1,0);
    // for(int i=0;i<1000;i++)
    for(int i=0;1;i++)
    {
        df = df*(a+double(i))*z/double((i+1)*(i+1));
        f = f + df;
        // if(abs(df.real())+abs(df.imag())<0.00001) break;
        if(abs(df.real())+abs(df.imag())<0.00000001) break;
    }
    return f;
}

///// 适用于|z|<45
complex<double> F_confluenthyper(complex<long double> a,complex<long double> z){
    complex<long double> f(1,0);
    complex<long double> df(1,0);
    // for(int i=0;i<1000;i++)
    for(int i=0;1;i++)
    {
        df = df*(a+(long double)(i))*z/(long double)((i+1)*(i+1));
        f = f + df;
        // if(abs(df.real())+abs(df.imag())<0.00001) break;
        if(abs(df.real())+abs(df.imag())<0.000001) break;
    }
    return complex<double>(f.real(),f.imag());
}




void MakeMomentumSpace_pLambda_V3(int seed);
void MakeMomentumSpace_pLambda_V4(int seed);
void MakeMomentumSpace_pLambda_V5(int seed);
void MakeMomentumSpace_pLambda_V6(int seed);
void MakeMomentumSpace_pLambda_V7(int seed);
void MakeMomentumSpace_pLambda_V2();
void MakeMomentumSpace_PP_V3();
void LL_pLambda_v0();
void LL_pLambda(int seed);
void LL_pLambda_All();
void LL_pLambda_All_seed(int seed);
void LL_pLambda_All_Multi_v0();
void LL_pLambda_All_Multi();

void LL_pp();

void LL_cpp(int seed=0){
    // MakeMomentumSpace_pLambda_V3(seed);
    // MakeMomentumSpace_pLambda_V4(seed);
    // MakeMomentumSpace_pLambda_V5(seed);
    // MakeMomentumSpace_pLambda_V6(seed);
    // MakeMomentumSpace_pLambda_V7(seed);//推荐
    // LL_pLambda_All();
    // LL_pLambda_All_Multi();
    LL_pLambda_All_seed(seed);
    // MakeMomentumSpace_pLambda_V2();
    // MakeMomentumSpace_PP_V3();
    // LL_pLambda();
    // LL_pp();

}

#include <thread>
using namespace std;

void LL_pLambda_All_Multi_v1(){

    ROOT::EnableThreadSafety();

    ///// scale5 <<<
	int nbin_f0=25;
    double f0_min = 1.5;
    double f0_max = 3.5;
	int nbin_d0=25;
    double d0_min = 0;
    double d0_max = 8;
	int nbin_R=25;
    double R_min = 1.8;
    double R_max = 3.8;
	int N_seed  = nbin_f0*nbin_d0*nbin_R;
    ///// scale5 >>>

    vector<thread*> thread_now;

    int N_Core = 10;

    for(int seed = 0;seed<N_seed;seed++){
        if(thread_now.size()<N_Core){
            auto this_thread = new thread(LL_pLambda_All_seed,seed);
            thread_now.push_back(this_thread);
        }
        if(thread_now.size()>=N_Core||seed==N_seed-1){
            // for(int i=0;i<thread_now.size();i++){
            for(auto it = thread_now.begin();it!=thread_now.end();it++){
                (*it)->join();
                delete *it;
                thread_now.erase(it);
            }
        }
    }

}


void LL_pLambda_v0(){

    // TChain* chain = new TChain("pLambda_Momentum");
    TChain* chain = new TChain("pLambda_Momentum_Iy0");
    // chain->Add("pLambda_Momentum.root"); // 600w
    // chain->Add("pLambda_Momentum_Iy0.small60w.root"); // 60w
    // chain->Add("pLambda_Momentum_Iy0.root"); // 6000w
    // chain->Add("pLambda_Momentum_Iy0.1wPair1Bin.root"); // 6000w
    chain->Add("pLambda_Momentum_Iy0.5kPair1Bin.3.root"); // 6000w
    // double *k = new double(2);
    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    
    // chain->SetBranchAddress("k", k);
    chain->SetBranchAddress("k", &k);
    // chain->SetBranchAddress("r", &beta);//beta
    chain->SetBranchAddress("beta", &beta);//beta
    // chain->SetBranchAddress("beta", &beta);//beta


    Int_t nevent = chain->GetEntries();
    cout<<Form("nevent=%d\n",nevent);


    TProfile* CF_LL;
    CF_LL = new TProfile("CF_pLambda", "CF_pLambda", 25, 0, 250);

    double f0 = 2.63;
    double d0 = 2.18;
    double R = 3.45;
    // f0=0.3;d0=0.5;R=2.3;
    double R_rel = R*sqrt(2);
    TLorentzVector r_rel;
    gRandom->SetSeed(4869);
    for(int n = 0;n<nevent;n++){
        chain->GetEvent(n);
        r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
        r_rel.Boost(*beta);
        double r_star = r_rel.Vect().Mag(); // fm
        double k_star = k->Mag()/hbar; // 1/fm

        double wei = 0;
        
        complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
        // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
        // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
        double kr = k_star*r_star;
        double k_r = r_rel.Vect()*(*k)/hbar;
        // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
        // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
        // wei = wei1+wei2;
        wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

        if(r_star<5*hbar){
            wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
        }

        CF_LL->Fill(k->Mag()*1000, wei);
    }
    TCanvas* C_CF = new TCanvas("C_CF","C_CF",1000,800);
    CF_LL->Draw();
    CF_LL->GetXaxis()->SetTitle("k/MeV");
    CF_LL->GetYaxis()->SetTitle("C(k)");
    CF_LL->SetMinimum(0.4);
    // CF_LL->SetMaximum(2);
    CF_LL->SetMaximum(3.6);
    CF_LL->SetAxisRange(0,200.5);
    // CF_LL->SetAxisRange(0,250.5);
    // CF_LL->SetAxisRange(0,0.25);
    C_CF->Print("CF_LL_pLambda.20230218.5k.png");
    TFile* CF_LL_File = new TFile("CF_LL_pLambda.20230218.5k.root","RECREATE");
    C_CF->Write();
    CF_LL->Write();
    CF_LL_File->Close();

}


#include <sys/time.h>
#include <stdio.h>
#include <unistd.h>

void LL_pLambda_All_seed(int seed){

    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    TTree*  MomentemTree[4];

    // auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");
    // auto file = new TFile("pLambda_Momentum_v7_5wPair1Bin.root","READ");
    auto file = new TFile("pLambda_Momentum_v6_5wPair1Bin.root","READ");

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        MomentemTree[I_y]->SetBranchAddress("k", &k);
        MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        // MomentemTree[I_y]->LoadBaskets(4000000000);
    }



    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale4")+TString(".root"),"RECREATE");
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale5")+TString(".root"),"RECREATE");

    

    // ///// scale3 <<<
    // int nbin_par = 50;
    // int N_seed = nbin_par*nbin_par*nbin_par;
    // double f0_min = 1.5;
    // double f0_max = 2.8;
    // double d0_min = 2.0;
    // // double d0_max = 7.0;
    // // double d0_max = 11.0;
    // double d0_max = 20.0;
    // double R_min = 1.0;
    // double R_max = 4.0;
    // ///// scale3 >>>

    // ///// scale4 <<<
	// int nbin_f0=20;
    // double f0_min = 0.3;
    // double f0_max = 4.3;
	// int nbin_d0=26;
    // double d0_min = 0.5;
    // double d0_max = 8.3;
	// int nbin_R=7;
    // double R_min = 2.3;
    // double R_max = 3.7;
	// int N_seed  = nbin_f0*nbin_d0*nbin_R;
    // ///// scale4 >>>

    ///// scale5 <<<
	int nbin_f0=25;
    double f0_min = 1.5;
    double f0_max = 3.5;
	int nbin_d0=25;
    double d0_min = 0;
    double d0_max = 8;
	int nbin_R=25;
    double R_min = 1.8;
    double R_max = 3.8;
	int N_seed  = nbin_f0*nbin_d0*nbin_R;
    ///// scale5 >>>

    // TProfile* CF_LL[N_seed][4];
    TProfile* CF_LL[4];
    // for(int seed = 0;seed<N_seed;seed++){
        cout<<"seed="<<seed<<endl;
        // struct timeval tv_begin, tv_end;
        // gettimeofday(&tv_begin, NULL);

        int f0_index = seed/(nbin_d0*nbin_R);
        int d0_index = (seed%(nbin_d0*nbin_R))/nbin_R;
        int R_index = (seed%(nbin_d0*nbin_R))%nbin_R;

        double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_f0);
        double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_d0);
        double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_R);


        

        TString Name = Form("CF_pLambda_f0_%.3f_d0_%.3f_R_%.3f",f0,d0,R);
        // TString CF_FileName = Form("pLambda_CF_v7_seed/All_CF_pLambda_v7_scale5_seed%d_",seed)+Name+TString(".root");
        // TString CF_FileName = Form("pLambda_CF_v7_seed/CF_pLambda_v7_scale5_seed%d_%s.root",seed,Name.Data());
        TString CF_FileName = Form("pLambda_CF_v6_seed/CF_pLambda_v6_scale5_seed%d_%s.root",seed,Name.Data());
        TFile* CF_LL_File = new TFile(CF_FileName,"RECREATE");
        cout<<CF_FileName<<endl;

        for(int I_y =0;I_y<4;I_y++){
            CF_LL[I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
        }


        Int_t nevent = MomentemTree[0]->GetEntries();
        cout<<Form("nevent=%d\n",nevent);



        
        // f0=0.3;d0=0.5;R=2.3;
        double R_rel = R*sqrt(2);
        TLorentzVector r_rel;
        gRandom->SetSeed(4869);
        for(int I_y =0;I_y<4;I_y++){
            for(int n = 0;n<nevent;n++){
                MomentemTree[I_y]->GetEvent(n);
                r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
                r_rel.Boost(*beta);
                double r_star = r_rel.Vect().Mag(); // fm
                double k_star = k->Mag()/hbar; // 1/fm

                double wei = 0;
                
                complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
                // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                double kr = k_star*r_star;
                double k_r = r_rel.Vect()*(*k)/hbar;
                // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei = wei1+wei2;
                wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

                if(r_star<5*hbar){
                    wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
                }

                CF_LL[I_y]->Fill(k->Mag()*1000, wei);
            }
        }

        // gettimeofday(&tv_end, NULL);
        // float process_interval   = (tv_end.tv_sec - tv_begin.tv_sec)*1000 + (float)(tv_end.tv_usec - tv_begin.tv_usec) /1000; //精度为毫秒
        // cout << "process_interval = "<<process_interval<<"ms"<<endl;
    // }

    // for(int seed = 0;seed<N_seed;seed++){
        for(int I_y =0;I_y<4;I_y++){
            CF_LL[I_y]->Write();
        }
    // }


    CF_LL_File->Close();

}


void LL_pLambda_All(){

    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    TTree*  MomentemTree[4];

    // auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");
    auto file = new TFile("pLambda_Momentum_v7_5wPair1Bin.root","READ");

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        MomentemTree[I_y]->SetBranchAddress("k", &k);
        MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        MomentemTree[I_y]->LoadBaskets(4000000000);
    }



    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale4")+TString(".root"),"RECREATE");
    TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale5")+TString(".root"),"RECREATE");

    

    // ///// scale3 <<<
    // int nbin_par = 50;
    // int N_seed = nbin_par*nbin_par*nbin_par;
    // double f0_min = 1.5;
    // double f0_max = 2.8;
    // double d0_min = 2.0;
    // // double d0_max = 7.0;
    // // double d0_max = 11.0;
    // double d0_max = 20.0;
    // double R_min = 1.0;
    // double R_max = 4.0;
    // ///// scale3 >>>

    // ///// scale4 <<<
	// int nbin_f0=20;
    // double f0_min = 0.3;
    // double f0_max = 4.3;
	// int nbin_d0=26;
    // double d0_min = 0.5;
    // double d0_max = 8.3;
	// int nbin_R=7;
    // double R_min = 2.3;
    // double R_max = 3.7;
	// int N_seed  = nbin_f0*nbin_d0*nbin_R;
    // ///// scale4 >>>

    ///// scale5 <<<
	int nbin_f0=25;
    double f0_min = 1.5;
    double f0_max = 3.5;
	int nbin_d0=25;
    double d0_min = 0;
    double d0_max = 8;
	int nbin_R=25;
    double R_min = 1.8;
    double R_max = 3.8;
	int N_seed  = nbin_f0*nbin_d0*nbin_R;
    ///// scale5 >>>

    TProfile* CF_LL[N_seed][4];
    for(int seed = 0;seed<N_seed;seed++){
        cout<<"seed="<<seed<<endl;
        struct timeval tv_begin, tv_end;
        gettimeofday(&tv_begin, NULL);

        int f0_index = seed/(nbin_d0*nbin_R);
        int d0_index = (seed%(nbin_d0*nbin_R))/nbin_R;
        int R_index = (seed%(nbin_d0*nbin_R))%nbin_R;

        double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_f0);
        double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_d0);
        double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_R);


        

        TString Name = Form("CF_pLambda_f0_%.3f_d0_%.3f_R_%.3f",f0,d0,R);

        for(int I_y =0;I_y<4;I_y++){
            CF_LL[seed][I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
        }


        Int_t nevent = MomentemTree[0]->GetEntries();
        cout<<Form("nevent=%d\n",nevent);



        
        // f0=0.3;d0=0.5;R=2.3;
        double R_rel = R*sqrt(2);
        TLorentzVector r_rel;
        gRandom->SetSeed(4869);
        for(int I_y =0;I_y<4;I_y++){
            for(int n = 0;n<nevent;n++){
                MomentemTree[I_y]->GetEvent(n);
                r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
                r_rel.Boost(*beta);
                double r_star = r_rel.Vect().Mag(); // fm
                double k_star = k->Mag()/hbar; // 1/fm

                double wei = 0;
                
                complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
                // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                double kr = k_star*r_star;
                double k_r = r_rel.Vect()*(*k)/hbar;
                // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei = wei1+wei2;
                wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

                if(r_star<5*hbar){
                    wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
                }

                CF_LL[seed][I_y]->Fill(k->Mag()*1000, wei);
            }
        }

        gettimeofday(&tv_end, NULL);
        float process_interval   = (tv_end.tv_sec - tv_begin.tv_sec)*1000 + (float)(tv_end.tv_usec - tv_begin.tv_usec) /1000; //精度为毫秒
        cout << "process_interval = "<<process_interval<<"ms"<<endl;
    }

    for(int seed = 0;seed<N_seed;seed++){
        for(int I_y =0;I_y<4;I_y++){
            CF_LL[seed][I_y]->Write();
        }
    }


    CF_LL_File->Close();

}



void LL_pLambda_All_Multi(){

    ROOT::EnableThreadSafety();

    // TVector3 *k = new TVector3(1,2,3);
    // TVector3 *beta = new TVector3(1,2,3);
    // TTree*  MomentemTree[4];
    TChain*  MomentemTree[4];

    // TString Momentum_filename("pLambda_Momentum_v6_5wPair1Bin.root");
    TString Momentum_filename("pLambda_Momentum_v7_5wPair1Bin.root");

    // auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");
    auto file = new TFile(Momentum_filename,"READ");

    for(int I_y =0;I_y<4;I_y++){
        // // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        // // MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        // MomentemTree[I_y] = (TChain*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        // // MomentemTree[I_y]->SetBranchAddress("k", &k);
        // // MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        // MomentemTree[I_y]->LoadBaskets(4000000000);


        MomentemTree[I_y] = new TChain(Form("pLambda_Momentum_Iy%d",I_y));
        // MomentemTree[I_y]->Add("pLambda_Momentum_v7/pLambda_Momentum.seed*.50Pair1Bin.V7.root");
        MomentemTree[I_y]->Add(Momentum_filename);
        // MomentemTree[I_y]->LoadBaskets(4000000000);
        // MomentemTree[I_y]->LoadBaskets(40000000000);
        
    }



    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale4")+TString(".root"),"RECREATE");
    TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale5")+TString(".root"),"RECREATE");
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/All_CF_pLambda_v6_scale5")+TString(".root"),"RECREATE");

    

    // ///// scale4 <<<
	// int nbin_f0=20;
    // double f0_min = 0.3;
    // double f0_max = 4.3;
	// int nbin_d0=26;
    // double d0_min = 0.5;
    // double d0_max = 8.3;
	// int nbin_R=7;
    // double R_min = 2.3;
    // double R_max = 3.7;
	// int N_seed  = nbin_f0*nbin_d0*nbin_R;
    // ///// scale4 >>>

    ///// scale5 <<<
	int nbin_f0=25;
    double f0_min = 1.5;
    double f0_max = 3.5;
	int nbin_d0=25;
    double d0_min = 0;
    double d0_max = 8;
	int nbin_R=25;
    double R_min = 1.8;
    double R_max = 3.8;
	int N_seed  = nbin_f0*nbin_d0*nbin_R;
    ///// scale5 >>>


    // const int N_Core = 40;//更慢
    // const int N_Core = 10;
    // const int N_Core = 20;
    const int N_Core = 15;
    // const int N_Core = 30;


    TProfile* CF_LL[N_seed][4];
    for(int seed = 0;seed<N_seed;seed++){
        cout<<"seed="<<seed<<endl;
        struct timeval tv_begin, tv_end;
        gettimeofday(&tv_begin, NULL);

        int f0_index = seed/(nbin_d0*nbin_R);
        int d0_index = (seed%(nbin_d0*nbin_R))/nbin_R;
        int R_index = (seed%(nbin_d0*nbin_R))%nbin_R;

        double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_f0);
        double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_d0);
        double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_R);


        

        TString Name = Form("CF_pLambda_f0_%.3f_d0_%.3f_R_%.3f",f0,d0,R);

        // for(int I_y =0;I_y<4;I_y++){
        //     CF_LL[seed][I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
        // }


        Int_t nevent = MomentemTree[0]->GetEntries();
        cout<<Form("nevent=%d\n",nevent);

        double R_rel = R*sqrt(2);

        for(int I_y =0;I_y<4;I_y++){
            auto workItem = [&](TTreeReader &reader) {
                TRandom3* thisRandom = new TRandom3(4869);
                // thisRandom->SetSeed(4869);
                TTreeReaderValue<TVector3> k(reader, "k");
                TTreeReaderValue<TVector3> beta(reader, "beta");
                // auto partialHisto = new TH1F("outHistoMP", "Random Numbers", 128, -4, 4);
                auto partialHisto = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
                while (reader.Next()) {
                    TLorentzVector r_rel;
                    r_rel.SetXYZT(thisRandom->Gaus(0,R_rel),thisRandom->Gaus(0,R_rel),thisRandom->Gaus(0,R_rel), 0);
                    r_rel.Boost(*beta);
                    double r_star = r_rel.Vect().Mag(); // fm
                    double k_star = k->Mag()/hbar; // 1/fm

                    double wei = 0;
                    
                    complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
                    // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                    // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                    double kr = k_star*r_star;
                    double k_r = r_rel.Vect()*(*k)/hbar;
                    // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                    // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                    // wei = wei1+wei2;
                    wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

                    if(r_star<5*hbar){
                        wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
                    }

                    partialHisto->Fill(k->Mag()*1000, wei);
                    // Info("WorkItem","wei=%f\n",wei);
                }
                return partialHisto;
            };
            ROOT::TTreeProcessorMP workers(N_Core);
            CF_LL[seed][I_y] = workers.Process(*(MomentemTree[I_y]), workItem, Form("pLambda_Momentum_Iy%d",I_y));
            CF_LL[seed][I_y]->Print();
        }
        
        gettimeofday(&tv_end, NULL);
        float process_interval   = (tv_end.tv_sec - tv_begin.tv_sec)*1000 + (float)(tv_end.tv_usec - tv_begin.tv_usec) /1000; //精度为毫秒
        cout << "process_interval = "<<process_interval<<"ms"<<endl;
    }

    for(int seed = 0;seed<N_seed;seed++){
        for(int I_y =0;I_y<4;I_y++){
            CF_LL[seed][I_y]->Write();
        }
    }


    CF_LL_File->Close();

}




void LL_pLambda_All_v1(){

    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    TTree*  MomentemTree[4];

    // auto file = new TFile("pLambda_Momentum.seed4.50Pair1Bin.V6.root","READ");
    // auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_5wPair1Bin.root","READ");
    // auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_v7_5wPair1Bin.root","READ");
    // auto file = new TFile("pLambda_Momentum_v7_2.5kPair1Bin.root","READ");
    // auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");
    auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        MomentemTree[I_y]->SetBranchAddress("k", &k);
        MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        MomentemTree[I_y]->LoadBaskets(4000000000);
    }



    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale2")+TString(".root"),"RECREATE");
    TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/All_CF_pLambda_v6_scale3")+TString(".root"),"RECREATE");

    int nbin_par = 50;
    int N_seed = nbin_par*nbin_par*nbin_par;

    ///// scale3 <<<
    double f0_min = 1.5;
    double f0_max = 2.8;
    double d0_min = 2.0;
    // double d0_max = 7.0;
    // double d0_max = 11.0;
    double d0_max = 20.0;
    double R_min = 1.0;
    double R_max = 4.0;
    ///// scale3 >>>

    // TProfile* CF_LL[8000][4];
    TProfile* CF_LL[N_seed][4];
    for(int seed = 0;seed<N_seed;seed++){
        cout<<"seed="<<seed<<endl;

        // double f0_min = 1.0;
        // double f0_max = 5.0;
        // double d0_min = 0.0;
        // double d0_max = 4.0;
        // double R_min = 1.0;
        // double R_max = 5.0;

        

        // int nbin_par = 20;
        // int nbin_par = 50;
        // int N_seed = nbin_par*nbin_par*nbin_par;

        ///// seed [0,8000) 20*20*20
        int f0_index = seed/(nbin_par*nbin_par);
        int d0_index = (seed%(nbin_par*nbin_par))/nbin_par;
        int R_index = (seed%(nbin_par*nbin_par))%nbin_par;

        double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_par);
        double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_par);
        double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_par);


        

        

        // double f0 = 2.63;
        // double d0 = 2.18;
        // double R = 3.45;

        TString Name = Form("CF_pLambda_f0_%.3f_d0_%.3f_R_%.3f",f0,d0,R);

        for(int I_y =0;I_y<4;I_y++){
            CF_LL[seed][I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
        }


        Int_t nevent = MomentemTree[0]->GetEntries();
        cout<<Form("nevent=%d\n",nevent);



        
        // f0=0.3;d0=0.5;R=2.3;
        double R_rel = R*sqrt(2);
        TLorentzVector r_rel;
        gRandom->SetSeed(4869);
        for(int I_y =0;I_y<4;I_y++){
            for(int n = 0;n<nevent;n++){
                MomentemTree[I_y]->GetEvent(n);
                r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
                r_rel.Boost(*beta);
                double r_star = r_rel.Vect().Mag(); // fm
                double k_star = k->Mag()/hbar; // 1/fm

                double wei = 0;
                
                complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
                // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                double kr = k_star*r_star;
                double k_r = r_rel.Vect()*(*k)/hbar;
                // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei = wei1+wei2;
                wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

                if(r_star<5*hbar){
                    wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
                }

                CF_LL[seed][I_y]->Fill(k->Mag()*1000, wei);
            }
        }

        // for(int I_y =0;I_y<4;I_y++){
        //     CF_LL[I_y]->Write();
        // }
    }

    for(int seed = 0;seed<N_seed;seed++){
        for(int I_y =0;I_y<4;I_y++){
            CF_LL[seed][I_y]->Write();
        }
    }
    // if(0){
    //     TCanvas* C_CF = new TCanvas("C_CF","C_CF",1000,800);
    //     CF_LL->Draw();
    //     CF_LL->GetXaxis()->SetTitle("k/MeV");
    //     CF_LL->GetYaxis()->SetTitle("C(k)");
    //     CF_LL->SetMinimum(0.4);
    //     // CF_LL->SetMaximum(2);
    //     CF_LL->SetMaximum(3.6);
    //     CF_LL->SetAxisRange(0,200.5);
    //     // CF_LL->SetAxisRange(0,250.5);
    //     // CF_LL->SetAxisRange(0,0.25);
    //     C_CF->Print("CF_LL_pLambda.20230218.5k.png");
    // }
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/")+Name+TString(".root"),"RECREATE");
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/")+Name+TString(".root"),"RECREATE");
    // C_CF->Write();

    // for(int I_y =0;I_y<4;I_y++){
    //     CF_LL[I_y]->Write();
    // }
    CF_LL_File->Close();

}


void LL_pLambda_All_v0(){

    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    TTree*  MomentemTree[4];

    // auto file = new TFile("pLambda_Momentum.seed4.50Pair1Bin.V6.root","READ");
    // auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_5wPair1Bin.root","READ");
    // auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_v7_5wPair1Bin.root","READ");
    // auto file = new TFile("pLambda_Momentum_v7_2.5kPair1Bin.root","READ");
    auto file = new TFile("pLambda_Momentum_v6_2.5kPair1Bin.root","READ");

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        MomentemTree[I_y]->SetBranchAddress("k", &k);
        MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        MomentemTree[I_y]->LoadBaskets(4000000000);
    }



    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/All_CF_pLambda_v7_scale2")+TString(".root"),"RECREATE");
    TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/All_CF_pLambda_v6_scale2")+TString(".root"),"RECREATE");

    for(int seed = 0;seed<8000;seed++){
        cout<<"seed="<<seed<<endl;

        // double f0_min = 1.0;
        // double f0_max = 5.0;
        // double d0_min = 0.0;
        // double d0_max = 4.0;
        // double R_min = 1.0;
        // double R_max = 5.0;

        double f0_min = 1.5;
        double f0_max = 2.8;
        double d0_min = 2.0;
        // double d0_max = 7.0;
        double d0_max = 11.0;
        double R_min = 1.0;
        double R_max = 4.0;

        int nbin_par = 20;

        ///// seed [0,8000) 20*20*20
        int f0_index = seed/(nbin_par*nbin_par);
        int d0_index = (seed%(nbin_par*nbin_par))/nbin_par;
        int R_index = (seed%(nbin_par*nbin_par))%nbin_par;

        double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_par);
        double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_par);
        double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_par);


        

        TProfile* CF_LL[4];

        // double f0 = 2.63;
        // double d0 = 2.18;
        // double R = 3.45;

        TString Name = Form("CF_pLambda_f0_%.3f_d0_%.3f_R_%.3f",f0,d0,R);

        for(int I_y =0;I_y<4;I_y++){
            CF_LL[I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
        }


        Int_t nevent = MomentemTree[0]->GetEntries();
        cout<<Form("nevent=%d\n",nevent);



        
        // f0=0.3;d0=0.5;R=2.3;
        double R_rel = R*sqrt(2);
        TLorentzVector r_rel;
        gRandom->SetSeed(4869);
        for(int I_y =0;I_y<4;I_y++){
            for(int n = 0;n<nevent;n++){
                MomentemTree[I_y]->GetEvent(n);
                r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
                r_rel.Boost(*beta);
                double r_star = r_rel.Vect().Mag(); // fm
                double k_star = k->Mag()/hbar; // 1/fm

                double wei = 0;
                
                complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
                // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
                double kr = k_star*r_star;
                double k_r = r_rel.Vect()*(*k)/hbar;
                // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
                // wei = wei1+wei2;
                wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

                if(r_star<5*hbar){
                    wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
                }

                CF_LL[I_y]->Fill(k->Mag()*1000, wei);
            }
        }

        for(int I_y =0;I_y<4;I_y++){
            CF_LL[I_y]->Write();
        }
    }
    // if(0){
    //     TCanvas* C_CF = new TCanvas("C_CF","C_CF",1000,800);
    //     CF_LL->Draw();
    //     CF_LL->GetXaxis()->SetTitle("k/MeV");
    //     CF_LL->GetYaxis()->SetTitle("C(k)");
    //     CF_LL->SetMinimum(0.4);
    //     // CF_LL->SetMaximum(2);
    //     CF_LL->SetMaximum(3.6);
    //     CF_LL->SetAxisRange(0,200.5);
    //     // CF_LL->SetAxisRange(0,250.5);
    //     // CF_LL->SetAxisRange(0,0.25);
    //     C_CF->Print("CF_LL_pLambda.20230218.5k.png");
    // }
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/")+Name+TString(".root"),"RECREATE");
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/")+Name+TString(".root"),"RECREATE");
    // C_CF->Write();

    // for(int I_y =0;I_y<4;I_y++){
    //     CF_LL[I_y]->Write();
    // }
    CF_LL_File->Close();

}

void LL_pLambda(int seed=0){

    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    TTree*  MomentemTree[4];

    // auto file = new TFile("pLambda_Momentum.seed4.50Pair1Bin.V6.root","READ");
    // auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_5wPair1Bin.root","READ");
    auto file = new TFile("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/pLambda_Momentum_v7_5wPair1Bin.root","READ");

    double f0_min = 1.0;
    double f0_max = 5.0;
    double d0_min = 0.0;
    double d0_max = 4.0;
    double R_min = 1.0;
    double R_max = 5.0;

    int nbin_par = 20;

    ///// seed [0,8000) 20*20*20
    int f0_index = seed/(nbin_par*nbin_par);
    int d0_index = (seed%(nbin_par*nbin_par))/nbin_par;
    int R_index = (seed%(nbin_par*nbin_par))%nbin_par;

    double f0 = f0_min + (f0_max - f0_min)*double(f0_index)/double(nbin_par);
    double d0 = d0_min + (d0_max - d0_min)*double(d0_index)/double(nbin_par);
    double R = R_min + (R_max - R_min)*double(R_index)/double(nbin_par);


    

    TProfile* CF_LL[4];

    // double f0 = 2.63;
    // double d0 = 2.18;
    // double R = 3.45;

    TString Name = Form("CF_pLambda_f0_%.2f_d0_%.2f_R_%.2f",f0,d0,R);

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = (TTree*)file->GetObjectChecked(Form("pLambda_Momentum_Iy%d",I_y),"TTree");
        MomentemTree[I_y]->SetBranchAddress("k", &k);
        MomentemTree[I_y]->SetBranchAddress("beta", &beta);
        CF_LL[I_y] = new TProfile(Name+Form("_Iy%d",I_y), Name+Form("_Iy%d",I_y), 25, 0, 250);
    }


    Int_t nevent = MomentemTree[0]->GetEntries();
    cout<<Form("nevent=%d\n",nevent);



    
    // f0=0.3;d0=0.5;R=2.3;
    double R_rel = R*sqrt(2);
    TLorentzVector r_rel;
    gRandom->SetSeed(4869);
    for(int I_y =0;I_y<4;I_y++){
        for(int n = 0;n<nevent;n++){
            MomentemTree[I_y]->GetEvent(n);
            r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
            r_rel.Boost(*beta);
            double r_star = r_rel.Vect().Mag(); // fm
            double k_star = k->Mag()/hbar; // 1/fm

            double wei = 0;
            
            complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
            // complex<double> c1 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
            // complex<double> c2 = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k*k,-k));
            double kr = k_star*r_star;
            double k_r = r_rel.Vect()*(*k)/hbar;
            // wei1 = 0.25*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
            // wei2 = 0.75*norm2(complex<double>(cos(kr),-sin(kr))+c1*complex<double>(cos(kr)/r,sin(kr)/r));
            // wei = wei1+wei2;
            wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );

            if(r_star<5*hbar){
                wei = norm2( complex<double>(cos(k_r),-sin(k_r)) );
            }

            CF_LL[I_y]->Fill(k->Mag()*1000, wei);
        }
    }
    // if(0){
    //     TCanvas* C_CF = new TCanvas("C_CF","C_CF",1000,800);
    //     CF_LL->Draw();
    //     CF_LL->GetXaxis()->SetTitle("k/MeV");
    //     CF_LL->GetYaxis()->SetTitle("C(k)");
    //     CF_LL->SetMinimum(0.4);
    //     // CF_LL->SetMaximum(2);
    //     CF_LL->SetMaximum(3.6);
    //     CF_LL->SetAxisRange(0,200.5);
    //     // CF_LL->SetAxisRange(0,250.5);
    //     // CF_LL->SetAxisRange(0,0.25);
    //     C_CF->Print("CF_LL_pLambda.20230218.5k.png");
    // }
    // TFile* CF_LL_File = new TFile(TString("pLambda_CF_v6/")+Name+TString(".root"),"RECREATE");
    TFile* CF_LL_File = new TFile(TString("pLambda_CF_v7/")+Name+TString(".root"),"RECREATE");
    // C_CF->Write();
    for(int I_y =0;I_y<4;I_y++){
        CF_LL[I_y]->Write();
    }
    CF_LL_File->Close();

}



void LL_pp(){
    // // MakeMomentumSpace(); // 慢，生成的文件大
    // // MakeMomentumSpace_V2(); // 推荐
    // MakeMomentumSpace_PP_V3(); // 推荐
    // return;

    TChain* chain = new TChain("pp_Momentum_Iy0");
    chain->Add("pp_Momentum_Iy0.5kPair1Bin.V3.root");
    TVector3 *k = new TVector3(1,2,3);
    TVector3 *beta = new TVector3(1,2,3);
    
    chain->SetBranchAddress("k", &k);
    chain->SetBranchAddress("beta", &beta);//beta


    Int_t nevent = chain->GetEntries();
    cout<<Form("nevent=%d\n",nevent);
    // nevent = 100000;

    TProfile* CF_LL;
    // CF_LL = new TProfile("CF_pp", "CF_pp", 25, 0, 250);
    CF_LL = new TProfile("CF_pp", "CF_pp", 80, 0, 200);
    TH1D* h_Wei;
    h_Wei = new TH1D("h_Wei","h_Wei",20,0,5);
    TH2D* h_k_Wei;
    // h_k_Wei = new TH2D("h_k_Wei","h_k_Wei",25,0,250,20,0,5);
    h_k_Wei = new TH2D("h_k_Wei","h_k_Wei",80,0,200,20,0,5);


    double f0_S0 = 7.8; // Singlet
    double d0_S0 = 2.77; // Singlet
    double f0_S1 = -5.4; // Triplet
    double d0_S1 = 1.7; // Triplet
    double R = 3;
    double R_rel = R*sqrt(2);
    TLorentzVector r_rel;
    gRandom->SetSeed(4869);

    // for(int n = 0;n<1;n++){
    for(int n = 0;n<nevent;n++){
        chain->GetEvent(n);
        r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);

        // {
        // // Test <<
        // // fsi_coor_:(x1,y1,z1,t1)=(-5.3116085348,0.7113185327,-1.8113493807,0.0000000000)
        // // fsi_coor_:(x2,y2,z2,t2)=(-0.4990404378,1.4457278382,-1.3983975247,0.0000000000)
        // // fsi_mom_:(px1,py1,pz1)=(1.0230957109,0.6744422544,-0.1692411230)
        // // fsi_mom_:(px2,py2,pz2)=(0.8710852864,1.0734581118,-0.1184039412)
        // TLorentzVector p1;
        // TLorentzVector p2;
        // p1.SetXYZM(1.0230957109,0.6744422544,-0.1692411230, 0.93827231);
        // p2.SetXYZM(0.8710852864,1.0734581118,-0.1184039412, 0.93827231);
        // r_rel.SetXYZT(-5.3116085348+0.4990404378,0.7113185327-1.4457278382,-1.8113493807+1.3983975247,0);

        // *beta = -(0.5*(p1+p2)).BoostVector();
        // // *beta = -(p1+p2).BoostVector();

        // p1.Boost(*beta);
        // p2.Boost(*beta);
        // *k = 0.5*(p1.Vect()-p2.Vect());
        // // Test >> wei = 1.2517640996397337 by Fortran | 1.251671 by C++
        // }

        r_rel.Boost(*beta);

        double r_star = r_rel.Vect().Mag(); // fm
        double k_star = k->Mag()/hbar; // 1/fm // k单位为GeV
        double kr = k_star*r_star;
        double k_r = r_rel.Vect()*(*k)/hbar;//k_dot_r
        if(abs(k_r)+abs(kr)>45) continue;// 超过这个范围合流超几何函数计算精度不够会溢出
        // double eta = alpha*Z*m/(hbar*(k->Mag())); // k,m,har均对应GeV //
        double eta = alpha*Z*m/(k->Mag()); // k,m,har均对应GeV //
        double Ac = TwoPi*eta/(exp(TwoPi*eta)-1);
        double h_eta = h_function(eta);
        double F,G;
        Coulomb_FG(0,eta,kr,F,G);

        complex<double> G_tilde = sqrt(Ac)*complex<double>(G,F);
        // complex<double> F_1 = F_confluenthyper(complex<double>(0,-eta),complex<double>(0,k_r+kr)); // \vec{r}
        // complex<double> F_2 = F_confluenthyper(complex<double>(0,-eta),complex<double>(0,-k_r+kr)); // -\vec{r}
        complex<double> F_1 = F_confluenthyper(complex<long double>(0,-eta),complex<long double>(0,k_r+kr)); // \vec{r}
        complex<double> F_2 = F_confluenthyper(complex<long double>(0,-eta),complex<long double>(0,-k_r+kr)); // -\vec{r}
        complex<double> PlaneWave_1 = complex<double>(cos(k_r),-sin(k_r));
        complex<double> PlaneWave_2 = complex<double>(cos(k_r),sin(k_r));
        complex<double> f_S0 = inverse(complex<double>(1/f0_S0+0.5*d0_S0*k_star*k_star-2*k_star*eta*h_eta,-k_star*Ac));
        complex<double> f_S1 = inverse(complex<double>(1/f0_S1+0.5*d0_S1*k_star*k_star-2*k_star*eta*h_eta,-k_star*Ac));
        double wei_S0 = 1./8.*Ac*norm2(PlaneWave_1*F_1+f_S0*G_tilde/r_star + PlaneWave_2*F_2+f_S0*G_tilde/r_star );
        double wei_S1 = 3./8.*Ac*norm2(PlaneWave_1*F_1+f_S1*G_tilde/r_star - (PlaneWave_2*F_2+f_S1*G_tilde/r_star) );

        // if(r_star<1){
        if(r_star<5*hbar){
            wei_S0 = 1./8.*Ac*norm2(PlaneWave_1*F_1 + PlaneWave_2*F_2 );
            wei_S1 = 3./8.*Ac*norm2(PlaneWave_1*F_1 - PlaneWave_2*F_2 );
        }

        double wei = wei_S0 + wei_S1;
        
        // if(wei>5||wei<0) continue;

        CF_LL->Fill(k->Mag()*1000, wei);
        h_k_Wei->Fill(k->Mag()*1000, wei);
    }
    TCanvas* C_CF = new TCanvas("C_CF","C_CF",1000,800);
    CF_LL->Draw();
    CF_LL->GetXaxis()->SetTitle("k/MeV");
    CF_LL->GetYaxis()->SetTitle("C(k)");
    CF_LL->SetMinimum(0);
    CF_LL->SetMaximum(2.5);
    CF_LL->SetAxisRange(0,200.5);
    C_CF->Print("CF_LL_pp.20230215.R3.5k.png");
    TFile* CF_LL_File = new TFile("CF_LL_pp.20230215.R3.5k.root","RECREATE");
    C_CF->Write();
    CF_LL->Write();
    h_k_Wei->Write();
    CF_LL_File->Close();

}




void MakeMomentumSpace_pLambda_V2(){

    TFile* MomentumTreeFile = new TFile("pLambda_Momentum_Iy2.5kPair1Bin.V2.root","RECREATE");
    TTree* MomentemTree = new TTree("pLambda_Momentum_Iy2","pLambda_Momentum tree");
    double k_star;
    // int I_y; // 0,1,2
    TVector3 beta;
    TVector3 k;
    MomentemTree->Branch("beta",&beta);
    MomentemTree->Branch("k",&k);

    auto CF_LL_count = new TProfile("CF_pLambda_count", "CF_pLambda_count", 25, 0, 0.25);
    int N_count4Bin = 5000;
    // int N_count4Bin = 10000;

    TLorentzVector p1;
    TLorentzVector p2;
    TLorentzVector p12;

    long int N = 1000000000000;//  
    gRandom->SetSeed(4869);
    for(long int n=0;n<N;n++){
    // for(double n=0;n<N;n++){
        if (n%1000000 ==0) cout<<"n="<<n<<",MomentemTree->GetEntries()="<<MomentemTree->GetEntries()<<endl;
        
        // if (int(n/10)%1000000 ==0) cout<<"n="<<n<<endl;

        double pt11 = gRandom->Uniform(0.5,2);
        double phi11 = TMath::TwoPi() * gRandom->Rndm();
        // double eta11 = gRandom->Uniform(-1,1);
        double eta11 = gRandom->Uniform(-1.85,0);
        p1.SetXYZM(pt11 * sin(phi11),pt11 * cos(phi11), pt11 * sinh(eta11), 0.938272);
        double y11 = p1.Rapidity() + 1.045;
        // if (y11 < 0 || y11 > 0.5) continue; // I_y0
        // if (y11 < 0.5 || y11 > 0.75) continue; // I_y1
        if (y11 < 0.75 || y11 > 1) continue; // I_y2

        double pt22 = gRandom->Uniform(0.5,2);
        double phi22 = TMath::TwoPi() * gRandom->Rndm();
        // double eta22 = gRandom->Uniform(-1,1);
        double eta22 = gRandom->Uniform(-1.85,0);
        p2.SetXYZM(pt22 * sin(phi22),pt22 * cos(phi22), pt22 * sinh(eta22), 1.11568);
        double y22 = p2.Rapidity() + 1.045;
        // if (y22 < 0 || y22 > 0.5) continue; // I_y0
        // if (y22 < 0.5 || y22 > 0.75) continue; // I_y1
        if (y22 < 0.75 || y22 > 1) continue; // I_y2


        // beta = -(0.5*(p1+p2)).BoostVector();
        beta = -(p1+p2).BoostVector();
        p1.Boost(beta);
        p2.Boost(beta);
        k = 0.5*(p1.Vect()-p2.Vect());
        k_star = k.Mag();
        if(k_star>0.25) continue;
        CF_LL_count->Fill(k_star,1);
        int count = CF_LL_count->GetBinEntries(CF_LL_count->FindBin(k_star));
        if(count>N_count4Bin) continue;

        MomentemTree->Fill();
        if(MomentemTree->GetEntries()>=N_count4Bin*25) break;
    }
    cout<<"MakeMomentumSpace Down,Get "<<MomentemTree->GetEntries()<<" Entries.\n";

    MomentumTreeFile->cd();
    MomentemTree->Write();
    MomentumTreeFile->Close();

    
}





void MakeMomentumSpace_pLambda_V6(int seed=0){
    // Add No y cut
    TFile fnew("3GeV_eta_pt.root");
    // TFile fnew("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/3GeV_eta_pt.root");
    // TFile fnew("/data4/qinz18/STAR/LL_Fit/LL_cpp/3GeV_eta_pt.root");
    TH1D *ptstr_p = (TH1D*)fnew.Get("hpt_p");
    TH1D *ptstr_la = (TH1D*)fnew.Get("hpt_la");

    double k_star;
    // int I_y; // 0,1,2
    TVector3 beta;
    TVector3 k;

    // TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum/pLambda_Momentum.seed%d.500Pair1Bin.V5.root",seed),"RECREATE");
    // TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum/pLambda_Momentum.seed%d.50Pair1Bin.V6.root",seed),"RECREATE");
    TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum_v6/pLambda_Momentum.seed%d.50Pair1Bin.V6.root",seed),"RECREATE");
    TFile* MomentumTreeFile[3];
    TTree* MomentemTree[4];
    TProfile* CF_LL_count[4];

    int N_Bin=25;
    int N_count4Bin = 50;
    // int N_count4Bin = 100;
    // int N_count4Bin = 500;
    // int N_count4Bin = 5000;
    // int N_count4Bin = 10000;

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = new TTree(Form("pLambda_Momentum_Iy%d",I_y),"pLambda_Momentum tree");
        MomentemTree[I_y]->Branch("beta",&beta);
        MomentemTree[I_y]->Branch("k",&k);
        CF_LL_count[I_y] = new TProfile(Form("CF_pLambda_count_Iy%d",I_y), Form("CF_pLambda_count_Iy%d",I_y), N_Bin, 0, 0.25);
    }

    

    TLorentzVector p1;
    TLorentzVector p2;
    TLorentzVector p12;

    long int N = 1000000000000;//  
    gRandom->SetSeed(4869+seed);
    for(long int n=0;n<N;n++){
    // for(double n=0;n<N;n++){
        if (n%1000000 ==0) cout<<"n="<<n<<",MomentemTree->GetEntries()="<<MomentemTree[0]->GetEntries()<<","<<MomentemTree[1]->GetEntries()<<","<<MomentemTree[2]->GetEntries()<<endl;
        
        // if (int(n/10)%1000000 ==0) cout<<"n="<<n<<endl;
        int I_y_p = -1;
        int I_y_lambda = -1;
        // double pt11 = gRandom->Uniform(0.5,2);
        double pt11 = ptstr_p->GetRandom();
        if(pt11 < 0.5 || pt11 > 2.0)continue;
        double phi11 = TMath::TwoPi() * gRandom->Rndm();
        // double eta11 = gRandom->Uniform(-1,1);
        double eta11 = gRandom->Uniform(-1.85,0);
        p1.SetXYZM(pt11 * sin(phi11),pt11 * cos(phi11), pt11 * sinh(eta11), 0.938272);
        double y11 = p1.Rapidity() + 1.045;
        // if (y11 < 0 || y11 > 0.5) continue; // I_y0
        // if (y11 < 0.5 || y11 > 0.75) continue; // I_y1
        // if (y11 < 0.75 || y11 > 1) continue; // I_y2
        if (y11 > 0 && y11 < 0.5) I_y_p=0;
        if (y11 > 0.5 && y11 < 0.75) I_y_p=1;
        if (y11 > 0.75 && y11 < 1) I_y_p=2;
        // if(I_y_p<0) continue;

        // double pt22 = gRandom->Uniform(0.5,2);
        double pt22 = ptstr_la->GetRandom();
        if(pt22 < 0.5 || pt22 > 2.0)continue;

        double phi22 = TMath::TwoPi() * gRandom->Rndm();
        // double eta22 = gRandom->Uniform(-1,1);
        double eta22 = gRandom->Uniform(-1.85,0);
        p2.SetXYZM(pt22 * sin(phi22),pt22 * cos(phi22), pt22 * sinh(eta22), 1.11568);
        double y22 = p2.Rapidity() + 1.045;
        // if (y22 < 0 || y22 > 0.5) continue; // I_y0
        // if (y22 < 0.5 || y22 > 0.75) continue; // I_y1
        // if (y22 < 0.75 || y22 > 1) continue; // I_y2
        if (y22 > 0 && y22 < 0.5) I_y_lambda=0;
        if (y22 > 0.5 && y22 < 0.75) I_y_lambda=1;
        if (y22 > 0.75 && y22 < 1) I_y_lambda=2;
        // if(I_y_lambda<0) continue;
        // if(I_y_lambda!=I_y_p) continue;
        // int I_y = I_y_lambda;
        int I_y=-1;
        if(I_y_lambda==I_y_p&&I_y_lambda>=0) I_y = I_y_lambda;



        // beta = -(0.5*(p1+p2)).BoostVector();
        beta = -(p1+p2).BoostVector();
        p1.Boost(beta);
        p2.Boost(beta);
        k = 0.5*(p1.Vect()-p2.Vect());
        k_star = k.Mag();
        if(k_star>0.25) continue;
        {
            int I_y=3;
            CF_LL_count[I_y]->Fill(k_star,1);
            int count = CF_LL_count[I_y]->GetBinEntries(CF_LL_count[I_y]->FindBin(k_star));
            if(count<=N_count4Bin){
                MomentemTree[I_y]->Fill();
            }
        }
        if(I_y>=0){
            CF_LL_count[I_y]->Fill(k_star,1);
            int count = CF_LL_count[I_y]->GetBinEntries(CF_LL_count[I_y]->FindBin(k_star));
            if(count>N_count4Bin) continue;
            MomentemTree[I_y]->Fill();
        }

        
        

        
        if((MomentemTree[0]->GetEntries()+MomentemTree[1]->GetEntries()+MomentemTree[2]->GetEntries()+MomentemTree[3]->GetEntries())>=4*N_count4Bin*N_Bin) break;
    }
    cout<<"MakeMomentumSpace Down,Get "<<MomentemTree[0]->GetEntries()<<"*3 Entries.\n";

    // for(int I_y =0;I_y<3;I_y++){
    //     MomentumTreeFile[I_y]->cd();
    //     MomentemTree[I_y]->Write();
    //     MomentumTreeFile[I_y]->Close();
    // }
    MomentumTreeFileAll->cd();
    MomentemTree[0]->Write();
    MomentemTree[1]->Write();
    MomentemTree[2]->Write();
    MomentemTree[3]->Write();
    MomentumTreeFileAll->Close();

    
}



void MakeMomentumSpace_pLambda_V7(int seed=0){ // CM
    // Add No y cut
    TFile fnew("3GeV_eta_pt.root");
    // TFile fnew("/star/u/qinz18/work/mywork/qinz18data/KF_LambdaProton/LLfit_ccode/3GeV_eta_pt.root");
    // TFile fnew("/data4/qinz18/STAR/LL_Fit/LL_cpp/3GeV_eta_pt.root");
    TH1D *ptstr_p = (TH1D*)fnew.Get("hpt_p");
    TH1D *ptstr_la = (TH1D*)fnew.Get("hpt_la");

    double k_star;
    // int I_y; // 0,1,2
    TVector3 beta;
    TVector3 k;

    // TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum/pLambda_Momentum.seed%d.500Pair1Bin.V5.root",seed),"RECREATE");
    // TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum/pLambda_Momentum.seed%d.50Pair1Bin.V7.root",seed),"RECREATE");
    TFile* MomentumTreeFileAll = new TFile(Form("./pLambda_Momentum_v7/pLambda_Momentum.seed%d.50Pair1Bin.V7.root",seed),"RECREATE");
    TFile* MomentumTreeFile[3];
    TTree* MomentemTree[4];
    TProfile* CF_LL_count[4];

    int N_Bin=25;
    int N_count4Bin = 50;
    // int N_count4Bin = 100;
    // int N_count4Bin = 500;
    // int N_count4Bin = 5000;
    // int N_count4Bin = 10000;

    for(int I_y =0;I_y<4;I_y++){
        // MomentumTreeFile[I_y] = new TFile(Form("pLambda_Momentum_Iy%d.seed%d.5kPair1Bin.V3.root",I_y,seed),"RECREATE");
        MomentemTree[I_y] = new TTree(Form("pLambda_Momentum_Iy%d",I_y),"pLambda_Momentum tree");
        MomentemTree[I_y]->Branch("beta",&beta);
        MomentemTree[I_y]->Branch("k",&k);
        CF_LL_count[I_y] = new TProfile(Form("CF_pLambda_count_Iy%d",I_y), Form("CF_pLambda_count_Iy%d",I_y), N_Bin, 0, 0.25);
    }

    

    TLorentzVector p1;
    TLorentzVector p2;
    TLorentzVector p12;

    TLorentzVector p_CM(0, 0, 1.1686, 1.498);
    TVector3 beta_CM = p_CM.BoostVector();

    long int N = 1000000000000;//  
    gRandom->SetSeed(4869+seed);
    for(long int n=0;n<N;n++){
    // for(double n=0;n<N;n++){
        if (n%1000000 ==0) cout<<"n="<<n<<",MomentemTree->GetEntries()="<<MomentemTree[0]->GetEntries()<<","<<MomentemTree[1]->GetEntries()<<","<<MomentemTree[2]->GetEntries()<<endl;
        
        // if (int(n/10)%1000000 ==0) cout<<"n="<<n<<endl;
        int I_y_p = -1;
        int I_y_lambda = -1;
        // double pt11 = gRandom->Uniform(0.5,2);
        double pt11 = ptstr_p->GetRandom();
        if(pt11 < 0.5 || pt11 > 2.0)continue;
        double phi11 = TMath::TwoPi() * gRandom->Rndm();
        // double eta11 = gRandom->Uniform(-1,1);
        double eta11 = gRandom->Uniform(-1.85,0);
        p1.SetXYZM(pt11 * sin(phi11),pt11 * cos(phi11), pt11 * sinh(eta11), 0.938272);
        double y11 = p1.Rapidity() + 1.045;

        
        p1.Boost(beta_CM);
        // cout<<"p1.Rapidity()="<<p1.Rapidity()<<","<<"y11="<<y11<<endl;
        
        // if (y11 < 0 || y11 > 0.5) continue; // I_y0
        // if (y11 < 0.5 || y11 > 0.75) continue; // I_y1
        // if (y11 < 0.75 || y11 > 1) continue; // I_y2
        if (y11 > 0 && y11 < 0.5) I_y_p=0;
        if (y11 > 0.5 && y11 < 0.75) I_y_p=1;
        if (y11 > 0.75 && y11 < 1) I_y_p=2;
        // if(I_y_p<0) continue;

        // double pt22 = gRandom->Uniform(0.5,2);
        double pt22 = ptstr_la->GetRandom();
        if(pt22 < 0.5 || pt22 > 2.0)continue;

        double phi22 = TMath::TwoPi() * gRandom->Rndm();
        // double eta22 = gRandom->Uniform(-1,1);
        double eta22 = gRandom->Uniform(-1.85,0);
        p2.SetXYZM(pt22 * sin(phi22),pt22 * cos(phi22), pt22 * sinh(eta22), 1.11568);
        double y22 = p2.Rapidity() + 1.045;

        p2.Boost(beta_CM);
        // cout<<"p2.Rapidity()="<<p2.Rapidity()<<","<<"y22="<<y22<<endl;

        // if (y22 < 0 || y22 > 0.5) continue; // I_y0
        // if (y22 < 0.5 || y22 > 0.75) continue; // I_y1
        // if (y22 < 0.75 || y22 > 1) continue; // I_y2
        if (y22 > 0 && y22 < 0.5) I_y_lambda=0;
        if (y22 > 0.5 && y22 < 0.75) I_y_lambda=1;
        if (y22 > 0.75 && y22 < 1) I_y_lambda=2;
        // if(I_y_lambda<0) continue;
        // if(I_y_lambda!=I_y_p) continue;
        // int I_y = I_y_lambda;
        int I_y=-1;
        if(I_y_lambda==I_y_p&&I_y_lambda>=0) I_y = I_y_lambda;



        // beta = -(0.5*(p1+p2)).BoostVector();
        beta = -(p1+p2).BoostVector();
        p1.Boost(beta);
        p2.Boost(beta);
        k = 0.5*(p1.Vect()-p2.Vect());
        k_star = k.Mag();
        if(k_star>0.25) continue;
        {
            int I_y=3;
            CF_LL_count[I_y]->Fill(k_star,1);
            int count = CF_LL_count[I_y]->GetBinEntries(CF_LL_count[I_y]->FindBin(k_star));
            if(count<=N_count4Bin){
                MomentemTree[I_y]->Fill();
            }
        }
        if(I_y>=0){
            CF_LL_count[I_y]->Fill(k_star,1);
            int count = CF_LL_count[I_y]->GetBinEntries(CF_LL_count[I_y]->FindBin(k_star));
            if(count>N_count4Bin) continue;
            MomentemTree[I_y]->Fill();
        }

        
        

        
        if((MomentemTree[0]->GetEntries()+MomentemTree[1]->GetEntries()+MomentemTree[2]->GetEntries()+MomentemTree[3]->GetEntries())>=4*N_count4Bin*N_Bin) break;
    }
    cout<<"MakeMomentumSpace Down,Get "<<MomentemTree[0]->GetEntries()<<"*3 Entries.\n";

    // for(int I_y =0;I_y<3;I_y++){
    //     MomentumTreeFile[I_y]->cd();
    //     MomentemTree[I_y]->Write();
    //     MomentumTreeFile[I_y]->Close();
    // }
    MomentumTreeFileAll->cd();
    MomentemTree[0]->Write();
    MomentemTree[1]->Write();
    MomentemTree[2]->Write();
    MomentemTree[3]->Write();
    MomentumTreeFileAll->Close();

    
}






void MakeMomentumSpace_PP_V3(){
    TFile* MomentumTreeFile = new TFile("pp_Momentum_Iy0.5kPair1Bin.V3.root","RECREATE");
    TTree* MomentemTree = new TTree("pp_Momentum_Iy0","pp_Momentum tree");

    TFile fnew("pt_0p4_2_proton_3GeV_Ur.root"); // 3 GeV
	TFile feta("eta_proton_3GeV_Ur.root"); // 3 GeV
	TH1D *ptstr = (TH1D*)fnew.Get("pt");
	TH1D *etastr = (TH1D*)feta.Get("eta");


    double k_star;
    // int I_y; // 0,1,2
    TVector3 beta;
    TVector3 k;
    // MomentemTree->Branch("k_star",&k_star);
    // MomentemTree->Branch("I_y",&I_y);
    // MomentemTree->Branch("r",&beta);
    MomentemTree->Branch("beta",&beta);
    MomentemTree->Branch("k",&k);

    int NBin = 80;
    double k_max = 0.2;//GeV

    // auto CF_LL_count = new TProfile("CF_pp_count", "CF_pp_count", 25, 0, 0.25);
    auto CF_LL_count = new TProfile("CF_pp_count", "CF_pp_count", NBin, 0, k_max);
    int N_count4Bin = 5000;
    // int N_count4Bin = 10000;

    TLorentzVector p1;
    TLorentzVector p2;
    TLorentzVector p12;

    // int N = 200000000;
    // int N = 1000000000;//Get 569700 Entries in 1000000000; 569700/1000000000. ~ 1/1755
    // double N = 10000000000;//  Down,Get 5695663 Entries.
    // double N = 100000000000;//  
    // long int N = 100000000000;//  
    long int N = 1000000000000;//  
    // TRandom3* random = new TRandom3(4869);
    gRandom->SetSeed(4869);
    for(long int n=0;n<N;n++){
    // for(double n=0;n<N;n++){
        if (n%1000000 ==0) cout<<"n="<<n<<",MomentemTree->GetEntries()="<<MomentemTree->GetEntries()<<endl;
        
        // if (int(n/10)%1000000 ==0) cout<<"n="<<n<<endl;

        double pt11 = ptstr->GetRandom();
        // double pt11 = gRandom->Uniform(0,2);
        double phi11 = TMath::TwoPi() * gRandom->Rndm();
        // double eta11 = gRandom->Uniform(-1,1); // -1.8-0?
        double eta11 = etastr->GetRandom(); //?
        p1.SetXYZM(pt11 * sin(phi11),pt11 * cos(phi11), pt11 * sinh(eta11), 0.938272);
        double y11 = p1.Rapidity() + 1.045;
        if (y11 < 0 || y11 > 0.5) continue; // I_y0
        // if (y11 < 0.5 || y11 > 0.75) continue; // I_y1
        // if (y11 < 0.75 || y11 > 1) continue; // I_y2

        // double pt22 = gRandom->Uniform(0,2);
        double pt22 = ptstr->GetRandom();
        double phi22 = TMath::TwoPi() * gRandom->Rndm();
        // double eta22 = gRandom->Uniform(-1,1);
        double eta22 = etastr->GetRandom();;
        // p2.SetXYZM(pt22 * sin(phi22),pt22 * cos(phi22), pt22 * sinh(eta22), 1.11568);
        p2.SetXYZM(pt22 * sin(phi22),pt22 * cos(phi22), pt22 * sinh(eta22), 0.938272);
        double y22 = p2.Rapidity() + 1.045;
        if (y22 < 0 || y22 > 0.5) continue; // I_y0
        // if (y22 < 0.5 || y22 > 0.75) continue; // I_y1
        // if (y22 < 0.75 || y22 > 1) continue; // I_y2


        // beta = -(0.5*(p1+p2)).BoostVector();
        beta = -(p1+p2).BoostVector();
        p1.Boost(beta);
        p2.Boost(beta);
        k = 0.5*(p1.Vect()-p2.Vect());
        k_star = k.Mag();
        if(k_star>k_max) continue;
        CF_LL_count->Fill(k_star,1);
        int count = CF_LL_count->GetBinEntries(CF_LL_count->FindBin(k_star));
        if(count>N_count4Bin) continue;

        MomentemTree->Fill();
        // if(MomentemTree->GetEntries()>=N_count4Bin*25) break;
        if(MomentemTree->GetEntries()>=N_count4Bin*NBin) break;
    }
    cout<<"MakeMomentumSpace Down,Get "<<MomentemTree->GetEntries()<<" Entries.\n";

    MomentumTreeFile->cd();
    MomentemTree->Write();
    MomentumTreeFile->Close();

    
}